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I.  INTRODUCTION 


This  study  is  the  first  step  in  determining  northing  values  in  the 
Universal  Transverse  Mercator  (UTM)  Military  Grid  System.  The  UTM  coordinate 
system  is  obtained  by  projecting' the  earth’s  surface  onto  an  oblate  cylinder 
slightly  smaller  than  the  earth's  radii.  The  northing  value  is  a  measure  of 
the  distance  from  the  equator  to  any  point  of  latitude  on  the  earth's  surface. 
The  first  step  in  determining  the  northing  value  requires  a  very  accurate 
calculation  of  the  distance  along  an  ellipse,  which  is  analogous  to  calcu¬ 
lating  the  UTM  northing  value  along  the  UTM  zone  central  meridian. 

The  purpose  of  this  study  was  to  determine  an  accurate  closed  loop  solu¬ 
tion  for  calculating  northing  values  along  the  central  meridian.  This  report 
derives  and  evaluates  three  different  techniques  for  determining  this  distance. 
They  are:  (1)  the  Gaussian  integration  of  a  second  order  elliptical  equation 
using  parametric  latitude,  (2)  the  closed  form  approximation  of  the  expanded 
and  integrated  elliptical  equation,  and  (3)  the  summation  of  average  radii  of 
curvature  as  a  function  of  the  geodetic  latitude  angle. 

A.  Derivation  of  the  Equation  for  an  Ellipse 

To  begin  this  study,  the  equation  of  an  ellipse  was  derived.  In  the 
following  drawing,  two  concentric  circles  of  radii  a,  b  with  centers  at  the 
origin  were  drawn.  Next,  a  ray  was  drawn  from  the  origin  cutting  both  circles 
at  the  points  Q  and  R  and  forming  the  angle  y.  P  is  then  the  intersection  of 
the  line  parallel  to  the  Y-axis  through  Q  and  the  line  parallel  to  the  X-axis 
through  R.  For  each  angle  y  there  is  determined  point  P,  whose  coordinates 
(X,  Y)  are  dependent  upon  y  . 


Since  00  =  a  and  OR  =  b,  X  and  Y  are  expressed  as  functions  of  v  by, 

X  =  a  cosy 
Y  =  b  sin  y 


1 


(1) 


From  the  equations  above,  it  can  be  seen  that 

X  Y 

—  =  cosy  and  —  =  siny 
a  b 


Then  by  squaring  and  adding  these  equations  together  and  using  the  trigono¬ 
metric  identity  cos^y  +  sin2y  =  l,  the  general  formula  for  an  ellipse  was 
derived.  Hence, 

x£  =  1 

a2  +  b2  (3) 

The  complete  ellipse  is  traced  out  as  y  (the  parametric  latitude)  goes 
from  0  to  2n.  Also  from  equation  (2)  and  the  same  trigonometric  identity  it 
was  found  that 

Y  b  _  (4) 

-  =  7  Tan  y 


Next,  by  looking  from  a  geocentric  perspective  the  tan  3  also  equalled 


Therefore  a  relationship  between  the  parametric  and  geocentric  angles  was 
established  by  combining  equations  (4)  and  (5). 


Tan  3  =  —  Tan  y  . 
a 


B.  The  Relationship  of  Beta  and  Gamma  to  Alpha 


Since  primary  emphasis  generated  around  geodetic  data,  a  relationship 
between  geodetic,  geocentric,  and  parametric  values  was  needed.  The  drawing 
below  (the  relationship  of  3  and  y  to  a)  shows  an  ellipse  and  the  geodetic, 
geocentric  and  parametric  angles. 


. X. 
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The  geodetic  angle,  a,  is  formed  when  taking  a  perpendicular  to  a  tangent 
at  any  point  on  Che  ellipse.  In  Equation  (3),  the  general  formula  for  an 
ellipse  was  differentiated  in  order  to  find  a  relationship  between  alpha  and 
the  beta  and  gamma  angles.  This  procedure  is  shown  below: 


Rearranging  Equation  (3)  into  differentiable  form  gives, 
X2b2  Y2a2  =  r 
a2b2  +  1^2  ~ 


y2b2  +  y2a2  =  a2b2  (constant). 
Differentiation  gives, 

2x  dx  b2  +  2y  dy  a2  =  0 


and 


2y  dy  a2  =  -2x  dx  b2 


(9) 


In  order  to  determine  the  tangent  of  alpha,  the  diagram  below  was  drawn. 


By  takmg  the  tangent  line  as  it  follows  the  arrow,  it  was  seen  that  dy  was 
positive  and  dx  was  negative,  thus  giving  a  negative  slope.  The  tangent  of 


alpha  was  found  by  rotating  a  90°  and  seeing  that  Tana 

Equation  (9)  for  — —  it  was  determined  that 
dy 


Tan  a  = 


-dx 

dy 


-dx 

dy 


So,  by  solving 


( in) 
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With  the  relationships  in  Equations  (6),  (9),  and  (10),  conversion 
equations  for  a  were  developed.  Substituting  the  equation  Tann  =  y/x  into 
Equation  (10)  gives 


and  by  substituting  Tanfl  =  b/ a  Tany  into  Equation  (10)  gives 


Equations  (11)  and  (12)  were  used  in  converting  geodetic  alpha  r  >  its 
geocentric  and  parametric  values. 

IT.  GAUSSIAN  INTEGRATION  OF  A  SECOND  ORDER  ELLIPTICAL  KOI'A  H' 'N 


Finally,  the  arc-length  formula  was  derived  from  tie*  ;  a r  in-- ’  :  i  **  ioi’i  ,os 
of  the  ellipse.  The  procedure  is  shown  as  follows: 


X  =  a  cosy 

Y  =  b  siny  ,  where  a  >  b. 


ds^  = 


dx2  +  dy2 


=  (a^sinZy  +  b2cos2y^dy2 
=  {^a^(l  -  co s 2 y^  +  b2cos2yjdy2 
r,  a2-b2  p  1 1/2  , 

Jds  =  aj  |\  -  —  yb—  cos2y|  c  dy 


ds  =  a 


i 


This  is  an  elliptical  integral  of  the  second  kind  where  e2 
the  eccentricity  of  the  ellipse.  Therefore, 


a2-b2 


Jds  =  aJ(l-e2cos2y^/2  dy 


,  where  e  is 


(14) 


The  first  of  the  three  methods,  the  Gaussian  numerical  integration,  using 
parametric  angles,  involved  using  the  integral  from  Equation  (14).  This 
integral  was  evaluated  from  lower  to  upper  gamma  in  steps  of  .5°  and  totaled 
for  a  cumulative  sum  of  arc-distance. 


III.  CLOSED  FORM  APPROXIMATION  OF  EXPANDED  AND  INTEGRATED  ELLIPTICAL 
EQUATION  -  METHOD  2 

Since  Equation  (14)  cannot  be  solved  in  closed  form,  the  binomial 
theorem  was  used  to  expand  the  integrand  into  an  integrable  closed  form  solu¬ 
tion.  By  using  the  definition  of  the  binomial  theorem 

(l+y)n  =  1+ny  +  n~~2 f ~ X  ■*■•••  etc.  (y2  <  i)f  (15) 
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and  setting 

o  o 

y  ~  e^cos^y 
n  =  V2, 

the  final  solution  became 

e2cos2  y  e^cos^  y 
1  2  "  8 


(16) 


(17) 


Therefore , 


/ds  -  - 


e2cos2y  e^cos^y  , 

—i  ~  dy  +  * 


and 


s  =  a 


’/ 

+  ei \ 

( 

1  4  64  J  Y  1 

{  8  32  J 

(18) 


(19) 


The  second  method  used  to  calculate  arc  distance  was  obtained  from 
Equation  (19)  which  was  evaluated  between  two  proper  limits.  Since  for  the 
case  under  study,  one  of  these  limits  is  y  =  0°,  Equation  (19)  can  be  evalu¬ 
ated  directly  as  a  function  of  y  (parametric  latitude). 


IV.  SUMMATION  OF  AVERAGE  RADIUS  OF  CURVATURE  -  METHOD  3 


The  third  and  last  method,  independently  derived  using  geodetic  lati¬ 
tudes,  involved  taking  the  radii  of  curvature  of  the  lower  and  upper  limits  (see 
drawing  below),  adding  them,  taking  their  average  and  multiplying  them  by  Act, 
where  Ax  is  the  change  in  geodetic  latitude. 


In  equation  form  it  is, 

s  ,  £  (ASkAJCU)  A«  , 

o 

where 

n  _  a( 1-e2) 

C  ”  (l-e2sin2'*)3/ 1 


(20) 


(21) 
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The  derivation  of  Equation  (21)  is  shown  in  the  Appendix. 


In  order  for  these  three  methods  to  be  calculated  smoothly  and  effi¬ 
ciently,  a  Fortran  program  called  ELIPSE  was  developed.  This  program,  as 
listed  below,  carries  out  all  the  aforenamed  integration  and  generates  a  table 
of  all  the  geodetic,  geocentric,  and  parametric  angles  along  with  the  three 
calculated  arc  distances  from  0°  to  90°  in  steps  of  .5°  (see  Table  1). 

This  study  was  the  first  step  in  determining  northing  values  along  the 
Central  Meridian.  Three  methods  analogous  to  chose  in  UTM  were  used  to  deter¬ 
mine  these  values.  The  first  method  was  a  numerical  integration  evaluation  of 
a  closed  form  integral.  The  second  method  was  a  series  approximation  of  the 
integral.  The  third  method  used  average  values  of  the  curvature  technique. 

The  results  from  comparisons  between  the  three  methods  show  differences 
of  Less  chan  a  meter.  The  conclusion  is  that  any  of  the  three  methods  could 
be  used  for  the  determination  of  UTM  northing  values  from  geographical  data 
without  accuracy  penalty. 


FORTRAN  PROGRAM 


C  PROGRAM  EL  IPSE 
C 

C  CALCULATES  THE  ADC LENGTH  BETWEEN  TWO  LIMITS  EXPRESSED  IN  ANGLES  P.Y  l!cE 
C  nr  the  GAUSSIAN  NUMERICAL  INTEGRATION  ROUTINE  .APPROXIMATE  CLOSED  FOL  •’ 

C  SOLUTION, AND  THE  AVERAGE  SUM  OF  THE  EQUATORIAL  RADII. 
r 

IMPLICIT  REAL*8  (A-H.0-7) 

EXTERNAL  EQUAT1 
C0W0N/AA/RADC.E2 
CHARACTERS  ANS 
C 

C  DATA  FOR  CALCULATIONS 

C 

DATA  A/63 78206. 400/ 

DATA  B/6356583. 8D0/ 

DATA  E2/.0067686580D0/ 

DATA  RADC/57. 20577951 3100/ 

DATA  PI2/1. 5707963267949000/ 

ADIF =0. 5/RADC 

EP2=( A**2-B**2)/B**2 

E4=E2**2 

FACI  =  ( B/A)  **2 

FAC2=B/A 

REWIND  3 

WRITE(2,100) 

100  FORMAT  (  '  ’  ,T7,  ’GEODETIC,  T21,  ’GEOCENTRIC ' ,  T38,  ’PARAMETRIC’ , 

•  Tfil, 'NUMERICAL ’,T75, ’APPROXIMATION' ,T°3, 'AVERAGE  RADII’ , Til 0, 
•'NITER' ,/,/,/) 

GO  TO  16 
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10  CONTINUE 
C 

C  DATA  CAN  BE  READ  IN  USING  F0R003.0AT  OR 

C  BY  TERMINAL  INPUT 

C 

C  WRITE! 6,*) 'DO  YOU  WISH  TO  COMTIMUE ?( Y /M) ' 

C  READ( 3, 1 01 )  ANS 

C  101  FORMAT! A) 

C  IF  (ANS.EQ. ’ Y ' )  THEN 

C  WRITE! 6,*) 'DO  YOU  WISH  TO  CONTINUE9! 1/2) ' 

READ! 3,*) IANS 
IF(IANS.EQ.l)  THEN 
GO  TO  16 
ELSE 

GO  TO  099 
EMD  IF 

16  CONTINUE 

C  WRITE! 6,*) ' INPUT  LOWER  LIMIT  VALUES  FOR  GEOOEDTIC  ALPHA  IN' 
C  WRITE! 6,*) 'DEGREES, MINUTES, AND  SECONDS. ’ 

READ! 3,*, END=550)AIDL, AIML, AISL 
C  WRITE! 6,*) 'POINTS9' 

READ! 3, *)M 

c  WRITE (6,*) 'SUB  INTERVALS?’ 

READ! 3 , *) MINT 

C  WRITE! 6,*) 'ACCURACY  FACTOR9' 

READ! 3, *)DEL 

C  WRITE!  G,*)  'MAX  INTEGRATION'S?  ' 

READ!  3,* MAX 


C 

C 

r 


c 

c 

r. 

r 

r 

C 

C 


r 

r 

c 


CHANGING  DEGREES, MINUTES, SECONDS  TO  DECIMAL  DEGREFS 
AND  CHANGING  DEGREES  To  RADIANS. 

CALL  DEGREE (DEGRE, AID L , AIML , AISL ) 

CALL  RADTANfRAD, DEGRE) 

ALOW=RAD 


INCREMENTING  UPPER  LIMIT  IN  STEPS  OF  .5  DEGREES 
AUP-ALOW+ADIF 

CORRESPONDING  GEOCENTRIC  AND  PARAMETRIC  LIMITS 
P(  ) GEOCENTRIC  G(  ) ^PARAMETRIC 

RLOW=OATAN( FAC1*DTAN( ALOW) ) 

GLOW=DATAM( FAC?*DTAN( ALOW) ) 

TEST  CASE  FOR  90  DEGREES 

IF  ( ALOW. GF. 1.55500)  THEN 
BUP=P  I? 

GUP = PI? 

ELSE 

B!JP  =  DATAM(  f  AC  1  *DTAN(  A;P ) ) 

GUP  =  DATAN!EAC?*-DTAN!A'!9)) 

ENO  IF 


C 


oooooo  oooo  o  o o 


C  GAUSSIAN1  INTEGRATION  .STEP  ONE. 

r 

CALL  GA!JSA(G, GDIF,  EQUAT1, GLOW  ,GUP,M,  MINT, DEL, MAX, MI  TER) 

TSUM=G*A 
SUM 1 =SUM 1 +TSUM 
C 

C  APPROXIMATE  CLOSEO  FORM  SOLUTION  .STEP  TVO. 

C 

EQUAT2=(  1 . 00- .  ?5D0*E2-3. DO/64  .  00*E4  )  *GUP-  ( 1  .r)0/8.D0*F?+1 . 00/3?.  00 
*  *E4  )  *DS  IN  (  2*GUP )  -  ( 1 . 00/256 . 90*E4  *DS  IN(  4*GU3 ) ) 

SUM2=EQ'JAT2*A 

C 

C  AVERAGE  SUM  EQUATORIAL  RAOII  .STE^  THREE. 

C 

RL=(A*( 1.00-E?) )/ 

■  (DSQRT(  1 .00- (E2*(0S IU(  ALO'.1) )  **2 ) )  *(  1 .  DO-  ( E2*(0SIM(  ALOV/) )  **.?))) 
RU=( A*( 1 . 00-E2 ) ) / 

■  (OSQRT(  1 . 00- ( E2*(0S II!(  AlIP ) )  VA2 ) )  *(  1 . DQ- ( E2 -(DSIM(  AUP ))**?)) ) 
DELTAA=AUP-ALOW 
EQUAT3=(  ( RL+P.ll )/? .  DO) *DELTAA 
SUM3=SUM3+EQUAT3 

RADIANS  BACK  TO  nEGREES, MINUTES, SECONDS 

CALL  TMVRAD(AID.AUP) 

CALL  INVRAD(BID.BUP) 

CALL  TMVRAD(GID.GUP) 

CALL  IMVDEG ( IAD,IAM,AS,AID) 

CALL  INVDEG( IBD, IBM,BS,BI0) 

CALL  INVDEG( IGD , IGM .GS.GID) 

OUTPUT 


WRITE(2, 200)  IAD.IAM.AS.IBD,  IBM, BS.IGO,  IGM  ,GS,SUM1,'SUM2,SUM3,  NITER 
200  FORMAT( '  ' ,T2, 13 ,T7 , 12, T10, F7 .4, T1 9, 13 , T24, 12, T27, F7. 4, T36, 13 ,T41 , 
■  I2.T44, F7.4, T53, FI  7. 0, T71 ,FI7. 0,T89,FI 7. 0,T107, 15,/ ,  / ) 

GO  TO  10 
550  CONTINUE 
909  CONTINUE 
STOP 
END 


SUBROUTINES  AND  FUNCTIONS 

DEGREE, RADI  AN, EQUAT1, INVDEG , INVRAD 

SUBROUTINE  DEGR EE ( OEGP. E , D EG , RMI N , SEC ) 

IMPLICIT  REAL*8  (A-H.O-Z) 
DEGRE=DABS(DEG)+RMIN/60.00+SEC/3600. DO 
DEGRE=DSIGN{DEGRE,DEG) 

RETURN 

END 


SUBROUTINE  RADIAM(RAf),DEGRE) 

COMMON/ AA/R  ADC,  E2 

RAD=DEGRE/RADC 

RETURN 

END 


FUNCTION  EQUAT1 (GAMMA) 

IMPLICIT  REAL*8  (A-H.O-Z) 

E2=. 0067686600 

EQUAT1 =D$ORT( 1 . 00- ( E2*0C0S(GAMMA) **2 ) ) 

RETURN 

END 


SUBROUTINE  INVRAD(DEGRE.RAD) 

IMPLICIT  REAL*8  (A-H.0-7) 

COMMON/ AA/R ADC, E2 

OEGRE=RAr'*RAQC 

RETURN 

END 


SUBROUTINE  INVOEG( I NT, MIN, SEC ,DEGRE) 
IMPLICIT  REAL*8  (A-H.O-Z) 

IMT= I ID INT ( DEGRE ) 

REAL=DFLOAT( INT) 

FRAC=DEGRE-REAL 

RSLT=FRAC*60.00 

M=IIDINT(P.SLT) 

REAL=DFLOAT(M) 

MIN=REAL 

FRAC=RSLT-REAL 

SEC=FRAC*60.D0 

RETURN 

END 


Angles  and  Calculated  Arc  Distances 
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APPENDIX 

From  a  standard  calculus  book*  Che  radius  of  curvature  is  defined  by  Che 
formula  below. 


[(f)2  +  (g’)2]3/2 
[f  g"  -  g*  f "I 


(  A-l ) 


Using  the  equations  from  the  ellipse  the  f  and  g  functions  are  defined. 


X  ■  f( y)  =  a  cos  y 
Y  =  g(  y)  =  b  sin  y 

where  y  is  the  parametric  latitude  defined  earlier. 


(A-2) 


Then  by  differentiating  these  equations  they  can  be  used  in  the  formula 
for  the  radius  of  curvature. 


f'  =  -a  sin  y 
g '  =  b  cos  y 


f”  =  -a  cos  y 
g"  =  -b  sin  y 

(a^sin^y  +  b2cos2y)2/2 
(ab  sin2y  +  ab  cos^y) 


=  — g-  [a2(l-cos2y)  +  b2cos2y)^2 
=  — r  [a2  -  (a2-b2)cos2 y] 2^2 


( A-3) 


RC  •  -ji  .2(1  -  cos2y  3/2 


-  a2-b2 

As  seen  earlier  in  this  report  e^  ■  — - —  .  By  transforming  this  equation 

a^ 

the  relationships  below  are  obtained. 


b^  2 

77  =  1  -  e2 


b  (l-e2) 


(A-4) 


0  3 

Therefore  ez  and  t-  can  be  substituted  into  the  Rq  equation. 


Rq  *  ( l-e2cos2 y)3/2 

=  a(l-e2cos2y)3/2 


*Calculus  with  Analytical  Geometry  by  Earl  W.  Swokowski 


RC  = 


a( 1-e^cos- y) 

aW'2 


(  A-5) 


Since  Che  radius  of  curvaturA  equation  is  needed  in  terms  of  a  (the  geode 
tic  latitude)  the  relationship  between  y  and  a  is  needed. 


Tany  =  —  Tan  a 
a 


If  a  triangle  like  the  one  shown  below  is  drawn,  then  the  cos^y  can  be 
obtained. 


y  (A-6) 

at  ■  i  -  i ...  —  -  .  . 

1  +  b^Tan^q  . 
a2 

Since  =  (l-e2)  then 


cos2  y 


1 _ 

l+( l-e2)Tan2q 


( A-7 ) 


By  expanding  Equation  (A-7)  and  substituting  it  into  Equation  (A-5)  then 
Equation  (21)  can  be  derived. 


cos2  y 


1 


1+Tan2  a_e2Tan2  q 


Rc 


a  1-e 


-p2 


3/2 


l+Tan2q-e2Tan2 i 


(  L-e2)  L/2 


_ Si _ _ \  3/2 

sec2 q-e^Tan/  ij 
(  l-e2 )  1/2 


i 

I 


a  1 


3/2 


1  l-e^sln^a 


cos^a 


( l-e^)  ^2 


l-e^sin2g-e^co6^a\  3/2 
l-e^sin^a 


(l-e2)V2 


a( l-e?) _ 

( l-e^sin2a)3/2 
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